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We systematically develop a density functional description for the equilibrium properties of a 
^ , two-dimensional, harmonically trapped, spin-polarized dipolar Fermi gas based on the Thomas- 

Q ■ Fermi von Weizsacker approximation. We pay particular attention to the construction of the two- 

dimensional kinetic energy functional, where corrections beyond the local density approximation 
must be motivated with care. We also present an intuitive derivation of the interaction energy 

■ functional associated with the dipolar interactions, and provide physical insight into why it can be 
C^") ■ represented as a local functional. Finally, a simple, and highly efficient self-consistent numerical 

procedure is developed to determine the equilibrium density of the system for a range of dipole 
^ . interaction strengths. 

a- 

6JT)f PACS numbers: 31.15.E-, 71.10.Ca, 03.75. Ss, 05.30.Fk 

03 " I. INTRODUCTION 

Si 

Ultra-cold, trapped dipolar quantum gases have received increasing attention over the past decade owing to the 
inherently interesting properties of the anisotropic, and long-range nature of the dipole-dipole interaction One of 
the important consequences of the anisotropy is that the interactions between the particles can be tuned from being 
I predominantly attractive to repulsive by simply changing the 3D trapping geometry, or for dipoles confined to the 2D 
" x-y plane, by adjusting the orientation of the dipoles relative to the z-axis [J Q • Therefore, novel physics in both the 
equilibrium and dynamic properties of such systems may be explored as a function of the strength of the interaction, 
the geometry of the confining potential, and the dimensionality of the system. 

While the degenerate dipolar Bose gas has been well studied experimentally and theoretically realizing a 
degenerate dipolar Fermi gas in the laboratory has proven to be much more elusive. One of the reasons for this 
is that the path to quantum degeneracy is impeded by the Pauli principle, which forbids s-wave collisions between 
identical atoms. Thus, early attempts to cool both magnetic and molecular dipolar Fermi gases below degeneracy were 
unsuccessful However, in the recent work of M. Lu et al. 0, this experimental hurdle was finally overcome, 

. resulting in the experimental realization of a spin polarized, degenerate dipolar Fermi gas. Specifically, using the 

■ method of sympathetic cooling, a mixture consisting of 161 Dy and the bosonic isotope 162 Dy, were cooled to T/Tp ~ 
0.2. In addition, this group was also able to evaporatively cool a single component gas of 161 Dy down to a temperature 
of T/Tp ~ 0.7. This latter result is presumed to arise from the rethermalization provided by the strong dipolar 

\ scattering between the 161 Dy atoms which have a large magnetic moment (/i ~ 10/is). 
7-H ■ The ability to fabricate such systems in the laboratory now opens the door for the investigation of both the 
equilibrium and dynamical properties of dipolar Fermi gases, and will enable contact to be made with the large body 
. i . of theoretical work already in the literature [l[. Moreover, it is now reasonable to expect that quasi-2D degenerate 
dipolar Fermi gases will also be realized experimentally, thereby allowing for studies into the stability, pairing, and 
superfluidity of low-dimensional dipolar systems, which to date have only been investigated theoretically 0, SQ. 

With a view towards ultimately calculating the collective mode frequencies, we develop in this paper a density- 
functional theory (DFT) for the equilibrium properties of a degenerate, harmonically trapped, spin polarized dipolar 
Fermi gas. Our theoretical framework is based on the Thomas-Fermi von Weizsacker (TFvW) approximation which 
was previously formulated in the context of degenerate electron gases (To| . The mathematical framework of the TFvW 
theory is very simple, numerically easy to implement, and computationally inexpensive. The TFvW theory has also 
been shown to provide an exceedingly accurate description of equilibrium properties, as well as collective excitations 
{i.e., magnetoplasmons) , of electronic systems in a variety of two-and three-dimensional confinement geometries (Tol — 
161 ] . Our purpose here is to take advantage of this approach, which is largely unknown in the cold-atoms community, 
and apply it to the dipolar Fermi gas. We will only address the equilibrium properties in this paper, and leave the 
presentation of the more involved mode calculation to a future publication. Moreover, in anticipation of forthcoming 
experiments, along with the goal of making contact with the recent theoretical work of Fang and Englert (l7j . we will 
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focus on the 2D geometry, although the extension of the theory to 3D is straightforward. The 2D geometry also allows 
for the development of exact analytical results, which we will exploit to test the efficacy of the TFvW approximation. 

The organization of our paper is as follows. In Sec. II, we construct the approximate kinetic energy functional 
for the trapped 2D system, and show that it is neccessary to go beyond the local-density approximation (LDA) in 
order to provide a more accurate determination of the ground state energy, along with physically reasonable density 
profiles, for the system. In Sec. Ill, we present an intuitive approach for the development of the interaction energy 
functional, in addition to providing physical insight into why the functional may be represented soley in terms of the 
local density. Section IV presents the details of the numerical procedure for implementing the TFvW, in addition to 
representative illustrations of the spatial density profile of the dipolar gas as the dipole-dipole interaction strength is 
changed. Finally, in Sec. V, we present our conclusions and closing remarks. 



II. KINETIC ENERGY FUNCTIONAL 

In the Kohn-Sham (KS) density-functional theory DFT [l8|, the noninteracting kinetic energy is treated exactly 
by solving ./V single-particle Schrddinger-likc equations, yielding the KS orbitals, from which the kinetic energy may 
be constructed. However, the KS DFT is not quite in keeping with the spirit of the Hohenberg-Kohn theorems, [19| 
which provide a mathematical justification for the solution to the many-body problem solely in terms of the density 
of the system (i.e., an orbital free formulation). Indeed, in its purest form, DFT has no need for the calculation 
of orbitals of any kind. However, an orbital free DFT requires the specification of a noninteracting kinetic energy 
density functional which is not known exactly. The purpose of this section is to provide an approximate, but accurate, 
kinetic energy density functional to be used in a DFT description of the ground state properties of a 2D harmonically 
trapped dipolar Fermi gas. 

The first level of approximation for the explicit construction of the kinetic energy density functional is the local- 
density approximation [20j . In this approach, the exact kinetic energy density for a homogeneous system is determined, 
after which the same expression is assumed to be true locally for the inhomogeneous system. The LDA is generally valid 
for systems that are only weakly inhomogeneous, but may be remarkably accurate even for strongly inhomogeneous 
systems [2l|, [22[ . In the case of a uniform 2D system, the noninteracting kinetic energy for a spin-polarized 2D Fermi 
gas is found to be 

k 

Here = 9(kp — k) is the zero-temperature Fermi occupation number, kp is the 2D Fermi wavevector given by 
kp = 4irfi2D, "2D is the uniform number density and A is the area of the system. Invoking the LDA, the 2D kinetic 
energy functional for an inhomogeneous system takes the form 

E K [n] = / dr— n(r) 2 . (2) 
J m 

If to this we add the energy Ep[n] of the particles interacting with an external potential w xt(r), we obtain the 
Thomas-Fermi (TF) energy functional 

dv^—n{Y) 2 + J drv cxt {v)n(r) . (3) 

A variational minimization of this equation with respect to the density leads to the Euler-Lagrange equation 

5Etf [n] 



6n(r) 



flTF = , (4) 



where the Lagrange multiplier htf (TF chemical potential) serves to fix the total number of particles. Using Eq. 
Eq. (Q| leads to the the TF spatial density, given explicitly as 

TTl 

n T F{r) = ^2 i^TF ~ WextOO) ■ (5) 



The density is seen to vanish on the surface defined by w ox t(r) = ^tf and is taken to be zero for all positions where 
^ext(r) > htf- This unphysical behaviour of the density is of course a consequence of the local nature of the TF 
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energy functional, and is also present in other spatial dimensions [23J. 

One may try to remove the unphysical behaviour of the TF density by improving upon the quality of the 2D kinetic 
energy functional. One possibility is introducing gradient corrections which take into account inhomogeneities of the 
density. However, when such a calculation is performed in 2D for a weakly inhomogeneous system, it is found that all 
gradient corrections vanish (l2l. l23l. I25U27I ] . The implication of this is that one cannot formally justify the inclusion of 
gradient corrections in 2D on the basis of a systematic expansion about the homogenous limit. This of course does 
not mean that the LDA term for the kinetic energy functional is exact since the TF approximation certainly does 
not generate the exact density of an inhomogeneous system. Nonlocal corrections to the kinetic energy functional are 
necessary, and it is plausible that they may still take the form of gradient corrections, similar to what are found in 
ID and 3D inhomogeneous systems. It is thus natural to ask if there is another way to at least motivate the inclusion 
of such gradient terms. 

In order to address this question, let us consider a 2D gas of noninteracting particles trapped within the harmonic 



confining potential v ox t(r) = mto^r 2 /2. It was shown by Brack and van Zyl 2l| that the exact spatial density is given 
by 



n cx (r) = - V(-1)"(M - n + l)L n (2r 2 )e~ r2 , (6) 

where L n (x) is a Laguerre polynomial [24j |. and all lengths are expressed in units of the harmonic oscillator length 
ai 10 = \J h/mujQ. The integer, M, counts the number of filled oscillator shells, i.e., the Fermi energy is given by 
Ep = hwo(M + 1). Integrating Eq. (|6]) over all space leads to the total number of particles in the system as a function 
of the shell index 

N(M) = \{M + 1)(M + 2) . (7) 



The exact total energy of this system is found to be 



hujQ 



E cx = -^NVT+8N . (8) 

By virtue of the equipartition of the kinetic and potential energies in a harmonic trap, the exact kinetic energy is 
given by 

6 
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N 3/2 + _J_ N l/2 + 0(^-1/2) 



(9) 



where in going to the last line in Eq. (l9l). we have made use of the large- N behaviour of the exact kinetic energy. 

The remarkable result found in Ref. [21| is that the LDA kinetic energy functional, Eq. (J2J) , integrates to the exact 
kinetic energy when the exact density is used as input. As mentioned above, this should not be misconstrued to mean 
that Eq. © is in fact the exact kinetic energy functional for the harmonically trapped system, i.e. that no corrections 
beyond the LDA are required. The exact spatial density n cx is emphatically not the density which variationally 
minimizes the TF energy functional. This density is given by Eq. ([5]) and may be written as 

«TF(r) = (R 2 TF - r 2 ) , (10) 



where Rtf = \J^^tf /mui 2 is the TF radius, and fiTF = \j2NtuxiQ. When Eq. (flQ|) is used in Eq. (O, the kinetic 
energy evaluates to 

EkVitf] = ^-hio^' 2 , (11) 

which is the leading term in the large- N expansion of ££. Thus, while Eq. © gives the quantum mechanical kinetic 
energy with the exact density as input, if the true variational density Eq. pU|) . is used instead, the resulting kinetic 
energy is always lower than the exact kinetic energy. The upshot of all of this is that the TF energy functional, 
Eq. ([3]), will always produce a kinetic energy that is lower than the true value. This implies that the kinetic energy 
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part of the TF functional has to be augmented by some correction in order to account for the second term in Eq. (|9|) . 

Motivated by the well-known nonlocal corrections found in ID and 3D, we take the correction to have the familiar 
von Wciszackcr (vW) form [23|, [28| 

E vW [n) = X vW (N)— / rfr |V } 1 , (12) 
8m J n(r) 

where \ v y/(N) is a parameter which in general can depend on the particle number, N. This functional has the 
following desirable properties: (i) it depends on the gradient of the spatial density and thus vanishes in the limit of 
a uniform system, (ii) it is positive definite, thus increasing the kinetic energy relative to the TF approximation, and 
(iii) it scales in the same way as EkIu] so that equipartition of kinetic and potential energy is preserved. The total 
energy functional for a noninteracting system in the TFvW approximation then reads (henceforth we suppress the 
iV-dependence of A v \y(-W)) 



E[n] = J dr 



nh 2 . , 2 x h 2 |Vn(r)| 2 

n r 2 + A vW ^ J 7~\ + «ext(r)n(r 

m 8m n(r) 



(13) 



The variational minimization of this energy functional is conveniently achieved by introducing the so-called von 
Weiszacker wave function ip(r) = ^Ju(y). The Euler-Lagrange equation then takes the form of a nonlinear Schrbdinger 
equation, 

h 2 

- A v w^VV(r) + v off (r)V(r) = fJ.tp(r) , (14) 
where the effective potential is given by 

fcft(r) = ^-V(r) 2 + Wext(r) . (15) 
m 

Since v e s (r) itself depends on ip(r), the solution of Eq. (|14[) must be determined self-consistently. The ground state 
solution ipo{ r ) with the normalization 

'dr|^ (r)| 2 = iV, (16) 

determines the self-consistent ground state density no(r) = V'o(r) 2 and the chemical potential [i. We now establish 
that the vW correction can account for the exact energy of the harmonically confined system with a parameter A v w 
which is weakly iV-dependcnt. 

For a given number of particles iV and a given value of A v \Vi we determine no(r) by solving the closed set of 
equations, viz., Eqs. ([T4]) and (fl"5j) . using the numerical method discussed in Sec. IV. This density is then used to 
evaluate E[hq] and the result is compared to E ex . The parameter A v w is then adjusted and the procedure is repeated 
until we achieve the equality 

E[n ] = E cx . (17) 

The net result of this procedure leads to the values of A v w plotted as a function of N in Fig. [TJ The ./V- dependence 
is indeed weak and suggests that A v w — 0.02-0.04 for N in the range 10 2 -10 6 . These values of A v w are considerably 
smaller than the value (~ 0.25) found to be appropriate in three-dimensions (29[. The inset to Fig. [T] illustrates the 
extrapolation to N — > oo and demonstrates that A v \y bas a nonzero limiting value. 

To sec in more detail how the vW energy accounts for the higher order terms in Eq. (j9)), it is convenient to expand 
E[no] in terms of the difference An = uq — utf- We have 

E[n ] = E TF [n ] + E vW [n ] 

= E TF [n T F\ + \mwl [ dr(r 2 - R^ F )An(r) + — [ rfr[An(r)] 2 + £ v w[n„] . (18) 

2 Jr>R TF m J 
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Equating (fT8|) to the large- TV expansion of Eq. (|SJ) and using Btf["tf] = {2^/2/3)huJoN 3 ^ 2 , we obtain 

lr>R TF 171 



^-^tHj N 1/2 ~ BvwM + ^mwg / dr{r 2 - R 2 TF )n (r) + f dr[n {r) - n TF (r)] 2 . (19) 



Each of the integrals on the right hand side, including the one defining £vw[no], have integrands that peak near 
r ~ Rtf- Thus the TV 1 / 2 term in Eq. which is a correction to the TF kinetic energy, can be interpreted as an 
edge correction. Since the last two terms on the right hand side of Eq. (|19|) are small in comparison to -E v w [no] , an 
approximate relation determining A v w would be 

E vW [no] = Avw^ / dr^^ ~ — L= h^N 1 ' 2 . (20) 
8m J n 12V2 

In applying this relation we again note that no in Eq. (|20p is an implicit function of A v w- Thus, the procedure 
described earlier is followed and A v w is adjusted until E v yv[no] is equal to the right-hand-side of Equation (|2"0"|) . This 
yields values of A v w which are about 10% larger than those obtained directly from Equation (fTTj) . 

To the extent that A v w is weakly TV-dependent, Eq. (f2"0")) indicates that the integral J dr|Vno| 2 /rio also scales roughly 
as N 1 / 2 . This iV-dependcnce is also exhibited by the integral with no replaced by n ex which is an indication that these 
densities are rather similar. In Fig. [2] we plot the self-consistent (solid line) and TFvW (dashed line) densities for 
different particle numbers. Even for the relatively small value of N = 231, these densities are in very good agreement 
and in fact are very close to Utf, except at the edge of the cloud. The differences between no and n ex are more clearly 
revealed by plotting i|Vn| 2 /n = \S7ip\ 2 in Fig.O The integrals of the curves shown in Fig. [3] typically differ by about 
15%. We note that the curves for the exact density (solid lines) exhibit prominent oscillations which are associated 
with the orbital shell structure. The fact that the TFvW curves (dashed lines) do not exhibit this shell structure is 
entirely expected in view of the semiclassical nature of the TFvW approximation 23] . However, we observe that the 
TFvW curves provide a smooth average of the shell oscillations (a well known feature of semiclassical theories (23j ) 
up to the edge of the cloud. At the edge, the TFvW approximation overestimates the exact value. In principle, this 
discrepancy can be reduced by including higher order gradient corrections to the kinetic energy functional, but for 
our purposes, this refinement is not neccessary. 

We summarize this section by noting that our analysis has shown that the corrections to the 2D LDA kinetic energy 
functional can indeed be represented in the gradient vW form, in spite of the fact that gradient expansion methods 
fail to produce any such terms in 2D. By comparing the results of the TFvW energy functional to the exact results 
for a harmonically confined noninteracting 2D gas, we are able to determine the vW parameter A v w and show that 
it is weakly A-dcpcndent. We cannot claim that these values of A v w are generally applicable for problems in 2D but 
we would expect them to be appropriate for densities which are similar to those of the harmonically confined system. 
In particular, we expect the vW functional to be applicable for a harmonically confined system in which interactions 
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FIG. 3: The quantity |V^| 2 = j\Wn\ 2 /n contributing to the von Weiszacker kinetic energy density in Eq. (fl2"l) . 
evaluated using the exact (solid line) and self-consistent (dashed line) spatial densities as a function of the radial 

distance r for different numbers of particles N . 



are also included, as discussed in the following section. 



III. DIPOLAR INTERACTIONS: HARTREE-FOCK APPROXIMATION 



Having developed the approximate kinetic energy functional in the previous section, we are now in a position to 
construct the energy functional which accounts for the dipolar interactions in a 2D spin polarized Fermi gas. The 
approach we adopt is essentially the one used for the analogous problem with Coulomb interactions in 2D degenerate 
electronic systems. However, as we shall see, dipolar interactions lead to a fundamentally different energy functional. 
At the level of the Hartree-Fock (HF) approximation [2(|, the energy of the spin polarized Fermi gas is Ehf = Ek+Em 
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where Ek = X)k n kfi k /2m is the noninteracting kinetic energy and the dipole interaction energy is 



E id = lj2 n * nk ' I( kk 'l^|kk') - (kk'|V dd |k'k) 



(21) 



kk' 



The first matrix element in Eq. pip is the direct term and the second is the exchange term. For the singular dipolar 
interaction of interest (i.e., with the spins polarized perpendicular to the plane), 



V dd {r) 



4 7rr 3 ' 



(22) 



each of these matrix elements are separately infinite. However, the sum of the two terms is finite as a result of the 
Pauli exclusion principle. This is seen most readily by writing the interaction energy as 



A f 

E dd = - J d 2 rV dd {v)n 2 2D g HF {r) 



where the HF radial distribution function is defined as 

nl D gHF(r) = n 2 



Evaluating Eq. 



we find 



9HF{r) = 1 - 



2J 1 (k F r) 
k F r 



-iW-r 



(23) 



(24) 



(25) 



where J\ (x) is the cylindrical Bessel function of order one |24| . The function gu F (r) — 1 has a range of approximately 
hp 1 and defines the 'exchange hole' commonly used in electron-gas theory. For small r, gHF(f) ~ j(k F r) 2 , and as a 
result, the integral in Eq. (f2"3"|) is well-behaved. We find 



1 



E dd = -Ahq/u, n 2D k F I dt— 



1 



2Ji(f) 
t 2 



A 64 2 _ 5/2 
A 45V^ 



(26) 



Equation (1261) has also been derived in the paper by Bruun and Taylor 0, and the final form is given by Fang and 
Englert [17] . In particular, Eq. (1261) suggests that one can define an interaction energy functional within a LDA having 
the following form: 



pLDA r i 

E dd [n\ 



d 2 r-C dd [niv)fl 2 : 
5 



(27) 



with C dd = (32/9v / 7r)/ioA i2 - In contrast to the Coulomb problem, this local functional presumably accounts for the 
total interaction energy. To see if this is indeed reasonable, one must investigate in more detail the effect of density 
inhomogeneities. 

Before doing so, we first provide an alternative derivation of Eq. (|26[) which assumes that V dd (r) has a well-defined 
Fourier transform (FT). This is achieved by defining a regularized dipole interaction V^ s (r) which does not have the 
r = singularity of V dd (r). As we shall see, a regularized interaction facilitates the corresponding analysis for an 
inhomogeneous system. Accepting for the moment that such an interaction is available, Eq. pip becomes 



K rcg 



d 2 k 



d 2 k' 



(2tt) 2 J (2tt) 



6(k F - k)8(k F - k') V£*(0) - V£ g (k - k') 



(28) 



where V^ s (k) is the two-dimensional FT of V dd s (r). With the change of variable k' = k q, Eq. p8p can be written 



^dd 



d 2 q 
(2tt)4 



^7(0) ~ ^f(q) d 2 k9(k F -k)e(k F -\k + q \)) 



(29) 
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The second integral is is just the area of overlap in momentum space of two circles of radius k F whose centres are 
separated by q. This area is given by 



A(q,k F ) = 2k 2 F 



cos 



Q 
2k F 



2k F 



q 

2k F 



9(2k F - q) 



(30) 



Substituting this result into Eq. (f2T)|). we find 



A k A 

reg _ Iih,p 
dd - 2tt 3 



dx: 



cos 1 x — xyl 



vr/(o)-vr/(2k FX ) 



(31) 



We now wish to show that this expression reduces to Eq. (|26[) using an appropriate limiting procedure. 

To this end, we must now specify the regularized dipole interaction. This is done by considering the interaction 
between two physical electric dipoles, each of which has a charge distribution of the form 



2pr r 2 
-e 



7T 3 / 2 CT 5 

The electrostatic interaction between two dipoles pi and P2 separated by r is 



(32) 



U(r) 



1 f rf3fc e ik-r(Pi__ k )(P2 ' k)e~ 



e J (2^) 3 

(Pi ■ V)(p 2 ■ V 
iiren 



k- 



lErf 
r 



V2a 



(33) 



where Erf (a;) is the error function. For r > <r, this interaction reduces to that of two point dipoles, varying as r~ 3 . 
However, for r <C cr, the interaction saturates at a constant value as a result of the overlap of the dipole charge 
distributions. 

Equation (|33[) can now be used to define a regularized magnetic dipole interaction for the spin-polarized Fermi gas 
by choosing Pi = P2 = pz, putting r = (x,y, 0) in Eq. (f3"3"]l . and replacing p 2 / £ o by /io/^ 2 - The regularized magnetic 
dipole interaction then reads 



v dd E ( x >y) = a*om 2 



47T 



d 3 k 
(2^)3 



r Erf 



J{k x x+k y y) k z -fc 2 cr 2 /2 



k 2 



V2cr 



2 1 



-r 2 /2a 2 



(34) 



Equation ([Mf approaches Eq. (|2^|) for r 3> a and saturates at no[i 2 /3(2tt) 3 > 2 ct 3 for r — > 0. The two-dimensional FT 
oiV* c /(x,y) is 



^(5) = MoM e 



dfc. 



27TCT 



2 2 

-9 o- /2 



27T fc 2 + g 2 
7T 



-fe 2 a 2 /2 



gcrErfc 



qa_ 

V2 



Vr/(0)~-K c d s (2k F x)^^k F x 



For ak F <C 1, we see that 



to leading order in ak F . Inserting this result into Eq. pip we obtain 

Hm = A 2^| 

o-^O dCl 45-7T 3 



(35) 



(36) 



(37) 



which is identical to Eq. (|26p . We thus see that the dipolar interaction energy can be obtained by taking the a — > 
limit in a calculation using an appropriately defined regularized dipole interaction. Of course the definition of the 
regularized potential is not unique, but the form we have chosen has particularly convenient properties. The essential 
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reason for being able to use this approach is that the final result is insensitive to the cutoff a when it becomes much 
smaller than the range kp 1 of the exchange hole. This calculation can be viewed as the momentum-space version of 
the real-space approach leading to Eq. (|23|) . 

We next proceed to a calculation of E dd for an arbitrary inhomogeneous system making use of the regularized dipole 
interaction defined above. As we shall see, this approach leads to a final result that is equivalent to that obtained in 
Reference [TtJ • The generalization of Eq. (|23[) to an inhomogeneous spin-polarized system is 

E M =\J d 2 r J d 2 r' [p(r, r)p(r', r') - p(r, r>(r', r)] V dd (r - r') , (38) 

where we have introduced the single particle density matrix defined as 

p(r,r') = ^^(r)^(r'). (39) 

2,OCC 

Here, the 4>i(r) are a set of single-particle states that correspond to a physical situation in which the density n(r) = 
p(r, r) is localized in space. The density matrix has the symmetry property p(r, r') = /o(r', r). The expression for E dd 
in Eq. (|38[) is well-defined even for the singular dipole interaction, however, it is more difficult to exhibit the explicit 
cancellation between the direct and exchange terms for an inhomogeneous system. To achieve this cancellation we 
make use of the regularized interaction in Eq. (|34j) which allows us to evaluate the the direct and exchange terms 
separately. The desired result is then obtained by taking the a — > limit at the end of the calculation. As we will 
show, the singular parts of the direct and exchange terms do in fact cancel exactly. 
The direct term is calculated most conveniently in momentum space. We have 



Efd' = \ J ^ 



To evaluate the exchange term, we introduce the centre-of-mass variable R = (r + r')/2 and the relative variable 
s = r — r'. The exchange term can then be written as 

Edl = -\j d 2 S VlT( S ) J d 2 R{p(R,s)} 2 , (41) 

where 

p(R,s) = p(R + ±s,R- |s) . (42) 
The symmetry property of p(r, r') implies that p(R, — s) = p(R, s). We now define the function 

/(b) = J d 2 R[p(R,s)] 2 , (43) 

which satisfies /(— s) = f(s), and write 

E ( dl - ~lJd 2 sVr/(s)f(s) 



From Eq. (|43p we have 



\ I d 2 sV d 7(s)[f(s) - f(Q)} \ I d 2 sVr/(s)f(0) . 



/(0) = I d 2 R[n(R)} 2 = J (0|n(q) 



(45) 



where the last equality follows from Parseval's theorem, and we have noted that n(R) = p(R, 0) = p(R, R). Combining 
Eq. g4]) and Eq. (gUJ), we obtain 

KTid) - ^7(0)] |n(q)| 2 - \ f d 2 s V^(s)[f( S ) /(0)] . (46) 



F _ 1 f d 2 q 

&dd — 



2 J (2tt) 2 
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It should be emphasized that this result is valid for any potential V dd s (r) which has a well-defined FT in the g->0 
limit. In case the potential does not have a well-defined q — > FT, one must revert to the expressions in Eqs. (|4"0")) 
and (|4Tj) . 

Making use of Eq. (pt5|) and taking the a — > limit, we obtain 



Edd = — 



= Ef}+E d \K (47) 

It can be shown that Eq. (|4"Tj) is identical to the result obtained by Fang and Englert [ItJ and for this reason we have 
adopted their notation for the two terms. The term E d ) appears explicitly in their paper but their E^J is given in a 
different form, being expressed in terms of the Wigner distribution function. 
The E dd term involves the quantity 

/(0)-/(s) = J d 2 R{[p(R,0] 2 -[p(R,s} 2 } 

d 2 i?[n(R)] 2 5 (s,R) . (48) 



Here we have defined the radial distribution function for the inhomogeneous system as 

^"I^F- < 49 » 

An approximation to E^J can be generated by making a local approximation for g(s,R), namely, 

g(s,Tl)~g HF (s;n{R)) , (50) 

where the radial distribution function of the uniform gas is evaluated for a density equal to n(R). With this replace- 
ment, we have 

E ( d l l - , (51) 

_____ 

as defined in Eq. (f2"T|) . The total interaction includes the manifestly nonlocal E dd term. 

We can check the validity of the LDA by evaluating Eq. (|47[) for a model inhomogeneous system. Specifically, 
we once again appeal to the ideal 2D Fermi gas confined by an isotropic harmonic trapping potential. This model 
is especially useful since the density matrix can be obtained in closed form for an arbitrary number of filled shells. 
Scaling all lengths by the harmonic oscillator length, a\ 10l the one-particle density matrix for M + 1 filled shells is 
given by [22j 



M 

p(R,s) = -J2(-V n Ln(2R 2 )e- R2 L\ I _ n ( S y2)e~ s2 / i , (52) 

7T * — ' 



n=0 



where L"(x) is the associated Laguerre polynomial [2J]. We observe that for this particular model system, p(R, s) 
depends only on the magnitudes of the vectors R and s, which is also true of all quantities derived from it. Putting 
s = in Eq. (|52p yields the density n(R) given by Eq. ([B]). (For convenience we drop the 'ex' subscript in the 
following.) The FT of n(R) is readily found to be given by 

n(q)=L 2 M (q 2 /2)e- q2/i ■ (53) 
We may also provide an explicit expression for the function /(s), namely, 

f(s) J d 2 R\p(R,s)\ 2 

M 

= ^EK^)] 2 ^ 2 ' ( 54 ) 

n=0 



11 



Since all of the functions required for the evaluation of Eq. (|47|) depend only on the magnitude of the coordinates, 
the angular integrations can be performed immediately, leading to 



^ 2,-^,2 MoM 2 [°°dsdf(s) 

ct/ *9 |n( ff )| / — 



*9 a l«(flr)l a -^S L / (55) 



*bo JO ™ho JO 

Using Eq. ([54]) along with Eq. (|53jl . we obtain 

,2 M 



p(i) MoM \ - 

dd " 87TO 3 ^ 



ho n=0 
2 roo 



ds 2^_ l(s2 / 2 )L r \( S 2 /2) e ^ /2 + / ds [L\^m 2 z- S l2 
o Jo 



(56) 



*ho JO 



Efl = -|gr / ^[^W^ 2 . (57) 



The function L 2 _ 1 (x), arising from the n = term in the summation of Eq. (|56[) . should be interpreted as zero. All of 
the integrals in Eqs. (f56|) and (|57l) can be evaluated analytically, using the following general result [22[ 



I mn (a,p,j) = / cte ^e^if^x)^) 
Jo 

T(l + a)r(n + 7 + l)r(/3 - a + m) 



F 2 (1 + a - p, -n, 1 + a; 1 + 7, 1 + a - P - m; 1) , (58) 



r(m + l)r(n + l)r(l + 7 )r(/3 - a) 
where 3-^2(0, c; <i, e; z) is the generalized hypergeometric function (24|. A direct application of Eq. (|58[) gives 

r 

"ho v 



F (l) _ MoM 1 \ - (n + l)r(n + 3/2) f4 p / 3 1 _. 9 1 _.,\ , P / 1 1 _. 9 1 „.,n /cqn 

A dd ~^3-y|L r(Tl+ 1) {3^ 3^2 ("2, 2'~ n ' 2 '-2 - "'Ij + 3^2 (-3, 2'~"' 2 '-2 - "'IJ) ' ( 5 9) 



.(2)_ MOM 2 1 r(M+|)(M + l)(M + 2) 3 

rfd ""^LiTl fpfTi) 3^(5=2.-^2.2-^.1) ■ (so) 

While these final expressions are not particularly illuminating, they do provide summations which are easily dealt 
with numerically. On the other hand, the integrands in Eqs. ()56j) and (|57j) become highly oscillatory when the index 
n is large and the evaluation of the integrals can be problematic without the use of specialized numerical integration 
techniques. 



N 


77.U) 771LDA 

^dd ^dd 


E (2) 
^dd 


AE/E% 


55 


54.4003 54.5724 


-6.8321 


15 


105 


168.6535 168.9366 


-15.3065 


10 


231 


670.199 670.718 


-40.9666 


7 


1326 


14266.3 14268.3 


-363.681 


3 


5151 


153345.6 153349.4 


-1983.06 


1 



TABLE I: Comparison of the dipolar interaction energies with the LDA. The last column corresponds to the relative 
percentage error between the exact total energy, E^ ct = Ew + E^J , and E^? A . The unit of energy is /io/i 2 / a ho- 

In Table U we give values of E^J , E^J and E^f A for a range of particle numbers, N. We see that E^J ~ E^ A to 
a very good approximation even for relatively small numbers of particles. To quantify this further we show in Fig. 0] 
a comparison of gHp(s;n(R)) and g(s,R) for M = 20 (N = 231) for R ranging from the centre of the cloud to its 
edge. We see that the local approximation in Eq. ([50)) is very good; the local Fermi wavevector fcp(i?) captures very 

nicely the extent of the exchange hole in the exact radial distribution function. These observations explain why E^J 
is so close to E d J d JA . Although these results have been obtained for the specific model of a harmonically confined gas, 
we would expect similar behaviour whenever the model density varies on a length scale which is large compared to 
the extent of the local exchange hole. 
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FIG. 4: Comparison of the exact radial distribution function (solid) with the local approximation (dashed). The 
values of R indicated are in units of the TF radius Rtf- The number of particles is N = 231 (M = 20). 



Table U also shows that the nonlocal contribution diminishes rapidly with respect to the local contribution with 
increasing N. Indeed, it is straightforward to show that the N ^s> 1 behaviour of the two contributions are E^) ~ TV 7 / 4 

and \E^J\ ~ iV 5 / 4 . Therefore, we find that \E^/E^\ ~ l/VN in the large- iV limit. As a result, the ioia/interaction 
energy is very well approximated by E^ A for large TV, as illustrated in Table |U This local energy functional for the 
total interaction energy can thus be trusted in applications, such as those typically encountered in traps, where the 
density of the system is a smooth and slowly varying function of position. 

We wish to emphasize that the locality of the interaction energy functional is a property of the dipolar interaction 
and is not generally valid. To illustrate this we can compare these results with those obtained for an interparticle 
Coulomb interaction, e 2 /Aire^r. There is no need to regularize the potential in this case and the direct and exchange 
terms can be evaluated directly from Equation (|38l) . For consistency we again consider a spin-polarized situation. 
The direct contribution is 

Eif = f- f 7^-Kq)l 2 , (6i) 

87re J (2n) z q 

where 2ir/q is the two-dimensional FT of r -1 . The exchange term reads 

Using Eqs. (|53|) and (|54j) for the isotropic two-dimensional harmonic trap, we find 

_e* 1 r(Af+f)r(A/ + 3) 

47re a h o 3\/2 T(M + l) 2 

and 



E& = ,_,„, - ^ V TV „ 2 /, \, 2 3 i ; a(-i,|,-Af;3,-f-M;l) , (63) 



47re a ho V2 ^ T(n + 1) 32[ 2 ' 2 ' '' 2 ' >' 1 } 

v n— 

The sum of these energies can be thought of as an approximation to the interaction energy of a 2D parabolic quantum 
dot. The 2D (spin-polarized) Dirac exchange functional in the LDA is given by 
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We compare the various energies in Table [TTJ We see that the LDA is again a very good approximation to the 
exchange energy. The direct Coulomb energy is of course inherently nonlocal and is seen to give the dominant 
interaction energy contribution for large N, in contrast to the situation for the dipolar interaction where the local 
contribution dominates. In fact, it can readily be shown that {E^/E^l ~ l/y/N for N > 1. 



N 




p (x),LDA 


E W 


AE/E% 


55 


-85.8544 


-85.4033 


683.615 


0.5 


105 


-192.347 


-191.746 


2119.36 


0.3 


231 


-514.802 


-513.958 


8421.97 


0.2 


1326 


-4570.15 


-4568.39 


179276.2 


0.04 


5151 


-24919.8 


-24911.7 


1.92699 xlO 6 


0.03 



TABLE II: Comparison of the Coulomb interaction energy 
LDA. The last column gives the relative percentage error 
exchange energy obtained from the Dirac functional, 



with the exchange energy calculated exactly and in the 
between the exact exchange energy, Eq. (|64p . and the 
Eq. ([6"5"]) . The energies are in units of e 2 /ineoaho- 



To summarize, we have shown that even for modest values of the particle number, N, the total interaction energy 
for the harmonically trapped, 2D dipolar Fermi gas can be approximately represented by a local energy functional. 
We must emphasize again that the locality of the density functional for the dipolar gas is a result of the short range 
~ 1/r 3 behaviour of the dipolc-dipole interaction, in contrast to e.g., the Coulomb interaction, in which case the total 
interaction energy cannot be reduced to a purely local form. 

We now have all of the necessary components to construct the TFvW energy functional for a 2D, harmonically 
trapped interacting dipolar Fermi gas, viz., 



E[n] = J dr 



l -C K n{vf + A vW ^ |V y 2 + lc dd [n{v)f' 2 + v cxt (r)n(r) 
I am n(r) 5 



(66) 



where Ck = 2nh 2 /m. Once again, introducing the vW wavefunction, and performing the variational minimization of 
Eq. (f66|) with respect to the density, gives 



- A vW — VV(r) + v eS (r)i>(r) = ntp(r) , 
2m 

where now, u e ff(r) also contains an interaction term, viz., 

v ce (r) = C K iP(r) 2 + C dd ^{rf + u ext (r) . 



(67) 



(68) 



Along with the normalization condition, Eq. (|16j) . Eq. (|67|) provides a complete description for the ground state spatial 
density of the system. 



IV. SELF-CONSISTENT EQUILIBRIUM SOLUTIONS 

The numerical self-consistent solution of Eq. (|67[) was achieved by means of imaginary time propagation (ITP) (30| . 
In this method, the time-dependent Schrodinger equation 

^>=-^W>, (69) 

is evolved in imaginary time, t — > — it, starting from some arbitrary initial state 0). The Hamiltonian governing 
the evolution is 

H(t) = -A vW ^ + v eS (r, r)=T + V(r) , (70) 
2m 
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where the time dependence arises from the dependence of v e g on the evolving density n(r, t) = V>(r, r) 2 . The evolution 
is carried out in a step-wise fashion according to 

iI>(t, t + At) = e- H ^ AT / h ij(r, r) . (71) 

The repeated application of the evolution operator yields a wave function ip{r,r) which eventually converges to the 
self-consistent ground state Vo(r) as r — > oo. 

The evolution in Eq. (|71|) is achieved by using the split-operator approximation [30l l3l| 

e -H(r)Ar/h _ g - V(r) At / 2h & -T At / h & -V (t) At / 2h ^ t^) 

together with fast Fourier transforms (FFT) [32j to convert between coordinate and momentum spaces. If a Cartesian 
grid is used in the x and y directions for our two-dimensional geometry, 2D FFTs are required. However, a more 
efficient algorithm is available if the system possesses circular symmetry. The kinetic energy operator in this case 
takes the form 

h 2 d 2 ft 2 ld_ 

T — -Avwt; rr + Avwt; -j- — T1+T2 . (73) 

2m ar z Ira r ar 

With this decomposition, the evolution operator becomes 

e -H( T )AT/H ^ e -V (t) At / 2h e -T 2 At / 2h er T! At /h e -T 2 At /2h e -V{ T ) At /2h /^N 

where the T\ step is again treated by means of FFTs, but now with respect to the one-dimensional radial variable r 
in the range [— R,R]. The kinetic step T2, on the other hand, is treated in coordinate space using a Crank- Nicholson 
method [321 ] which requires the solution of a tridiagonal system, viz., 

- <f>(r i+ i) + ai<j>(ri) + 4>{ri-i) = 4>{r i+1 ) + ai<f>(n) - , (75) 

4t*2 Ar - 

where a, = — , and 0(r,) is the wave function prior to the application of the Ti evolution operators. The 

A v w_At 

solution to Eq. (|75[) is uniquely determined by Dirichlct boundary conditions at the ends of the r, mesh, where the 
wave function is required to vanish. At the end of each time step we update H with the new n(r, r) which is properly 
normalized to N. The convergence criterion for achieving self-consistency is taken to be iV'fai r n) — "0( r ij T n-i)l < £ i 
where typically e < 10 -6 . Once self-consistency has been achieved, the chemical potential is given by [1 = (ipo\H \tpo) 
and the ground state energy is obtained from Equation (|66[) . This numerical procedure can also be adapted with 
minor modifications to spherically symmetric 3D systems. 
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FIG. 6: Density distributions for N = 1000. The curves with increasing radial extent correspond to Cdd/Cx = 0, 
0.2, 0.5, 1.0, 2.0, 5.0, 10.0. The inset shows the effective potential in the case of strong interactions, Cdd/Cx — 10.0. 

In Fig. [5] we present the TFvW self-consistent ground state density profile for N = 100 for various strengths of the 
dipolc-dipole interaction. With increasing interaction strength, the radius of the atomic cloud increases and the central 
density decreases, which is expected given the repulsive nature of the dipolc-dipole interaction for the spin-polarized 
case. As found for the Cdd = case considered earlier, the vW gradient term has the effect of giving a density that 
decays smoothly in the classically forbidden region defined by v c s(r) > [i. We recall that if the vW term is absent, 
the spatial density has an unphysical sharp cut-off at a radius R defined by v e s(R) = //, where // is the chemical 
potential in this approximation. For Cdd/Cx ~> 1, Roc TV 3 / 10 . In Fig. [6l we show the spatial density for N = 1000. It 
is apparent that the effect of the vW term becomes less significant as the number of particles is increased. This trend 
increases with increasing N and the density approaches the distribution found in the local density approximation. 
In the inset to both Figs. [5] and [5] we have also included representative plots of the effective potential for the case 
Cdd/CK = in Fig. [5] and Cdd/Cx = 10 in Fig. [6j The potential is essentially flat up to the edge of the cloud where 
it falls below the chemical potential and then rises quadratically. The main point to be taken from these curves is 
that the shape of the effective potential is not very sensitive to the introduction of interactions. Finally, we note that 
our equilibrium density profiles are similar to those found in Ref. \17\ . with the primary difference being that our 
self-consistent solutions exhibit a smooth decay into the classically forbidden region as a result of our inclusion of the 
vW gradient correction to the 2D Thomas-Fermi kinetic energy functional. 

V. CONCLUSIONS 

We have presented a mathematically simple, and computationally efficient DFT formulation of the equilibrium 
properties of a 2D trapped dipolar Fermi gas based on the Thomas-Fermi von Weizsacker approximation. One of 
the key elements of this work is the development of a kinetic energy functional appropriate for an inhomogeneous 
2D system. Specifically, we have shown that the addition of a vW-like gradient correction to the TF kinetic energy 
functional is needed in order to accurately determine the ground state energy, as well as to provide a physically 
reasonable density distribution at the edge of the cloud. While conventional gradient expansions in 2D fail to yield 
gradient corrections to the TF kinetic energy functional, our detailed analysis has provided a compelling argument 
for the inclusion of a vW-like gradient correction. We have also provided an alternative derivation of the interaction 
energy functional first obtained by Fang and Englert [13] • We find in particular that the exchange hole is a useful 
concept in formulating the local density approximation, and furthermore explains the underlying reasons behind the 
local nature of the total interaction functional. 

We have also presented a highly efficient self-consistent numerical scheme for determining the equilibrium spatial 
density distributions within the TFvW formalism. Our calculations illustrate how the strength of the (repulsive) 
dipolc-dipole interaction affects the size of the cloud and its spatial distribution. Although the vW gradient correction 
does not modify substantially the density in the interior of the cloud, it does yield a density which decays smoothly 
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into the classically forbidden region at the edge. This feature is particularly important in performing calculations of 
the collective mode frequencies using TFvW hydrodynamics [ToL fl6| . and will be addressed in a future publication. In 
addition, it will be of interest to investigate the affect of the anisotropy of the dipolc-dipole interaction in a confined 
2D system by considering the situation in which the spins are canted at some angle with respect to the 2D plane. 
Finally, although we have focused on the 2D geometry, we wish to emphasize that the extension of the TFvW theory 
to 3D is straightforward. 
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